GSEA analysis of SARC associated genes

Data Downloading

sratools

We downloaded the data using SRAtools on grace using the scripts generated by the following codes. We downloaded the sratools from github and installed it on grace with setup. All files downloaded were not zipped so we also zipped each file into fastq.gz file.

id.list<-readLines("/Users/yanxiting/Downloads/SRR_Acc_List(2).txt")
temp1<-paste("/home/xy48/scratch_palmer/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin/prefetch.3.0.0 -X 9999999999999 ",id.list,"\n",sep="")
temp2<-paste("/home/xy48/scratch_palmer/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin/fastq-dump.3.0.0 --split-files ./",id.list,".sra\n",sep="")
temp3<-rep("cd /home/xy48/scratch_palmer/SARC_10x/test/sra\n",length(temp1))
temp4<-paste("gzip ./",id.list,"_*.fastq\n",sep="")
temp5<-rep("ssh transfer\n",length(temp1))
cmd.out<-cbind(temp5,temp1,temp3,temp2,temp4)
temp.cmd<-apply(cmd.out,1,paste,collapse="",sep="")
cat(temp.cmd,file="./download_commands.txt",append=F,sep="")

changing names

The cell ranger pipeline recognize fastq file names are formated as [Sample Name]S1_L00[Lane Number][Read Type]_001.fastq.gz. There are 3 files per sample representing I1, R1 and R2 downloaded by sratools, which were named as *_1, *_2 and *_3.fastq files. To correctly run cell ranger on these files, we need to change the fastq.gz names.

# this was run on grace.
source("~/Rprogram/my_functions.R")
data.dir<-"/home/xy48/scratch_palmer/SARC_10x/test/sra"
filenames<-list.files(data.dir)
temp<-sapply(filenames,my.element.extract,splitchar="\\.",index=3)
filenames<-filenames[temp=="gz"]
filenames<-filenames[!is.na(filenames)]

for(i in 1:length(filenames)){
  
  if(substr(filenames[i],1,4)=="SRR9"){
    from.filename<-filenames[i]
    temp<-my.element.extract(filenames[i],splitchar="\\.",index=1)
    temp<-my.element.extract(temp,splitchar="_",index=2)
    if(temp=="1"){
        to.filename<-paste0(my.element.extract(filenames[i],splitchar="_",index=1),"_S1_L001_I1_001.fastq.gz")
    }
    if(temp=="2"){
        to.filename<-paste0(my.element.extract(filenames[i],splitchar="_",index=1),"_S1_L001_R1_001.fastq.gz")
    }
    if(temp=="3"){
        to.filename<-paste0(my.element.extract(filenames[i],splitchar="_",index=1),"_S1_L001_R2_001.fastq.gz")
    }
  
    file.rename(from=file.path(data.dir,from.filename),to=file.path(data.dir,to.filename))
    cat("from=",file.path(data.dir,from.filename),"\n","to=",file.path(data.dir,to.filename),"\n",sep="")
  }else{
    
    from.filename<-filenames[i]
    temp<-my.element.extract(filenames[i],splitchar="\\.",index=1)
    temp<-my.element.extract(temp,splitchar="_",index=2)
    if(temp=="1"){
        to.filename<-paste0(my.element.extract(filenames[i],splitchar="_",index=1),"_S1_L001_R1_001.fastq.gz")
    }
    if(temp=="2"){
        to.filename<-paste0(my.element.extract(filenames[i],splitchar="_",index=1),"_S1_L001_R2_001.fastq.gz")
    }
    if(temp=="3"){
        to.filename<-paste0(my.element.extract(filenames[i],splitchar="_",index=1),"_S1_L001_I1_001.fastq.gz")
    }
  
    file.rename(from=file.path(data.dir,from.filename),to=file.path(data.dir,to.filename))
    cat("from=",file.path(data.dir,from.filename),"\n","to=",file.path(data.dir,to.filename),"\n",sep="")    
    
  }
}

cell ranger

We rsync the fastq.gz files onto farnam under ~/scratch60/SARC_10X/fastq and run cell ranger on the files using commands generated by the following codes.

id.list<-readLines("/Users/yanxiting/Downloads/SRR_Acc_List(2).txt")
temp1<-paste("#!/bin/bash\n#SBATCH --time=24:00:00  --ntasks=1 --partition=general --cpus-per-task=8 --mem=80GB --job-name=",id.list," -o /home/xy48/scratch60/SARC_10X/cellranger_scripts/",id.list,".sh.o%J -e /home/xy48/scratch60/SARC_10X/cellranger_scripts/",id.list,".sh.e%J\n",sep="")
temp2<-rep("module load cellranger/5.0.0\nmodule load bcl2fastq2/2-20-0-foss-2018b\n",length(id.list))
temp3<-paste("cd /home/xy48/scratch60/SARC_10X/cellranger_results\n")
temp4<-paste("cellranger count --id=",id.list," --fastqs=/home/xy48/scratch60/SARC_10X/fastq --transcriptome=/home/xy48/scratch60/SARC_10X/refdata-gex-GRCh38-2020-A --localcores=8 --localmem=60 --sample=",id.list,"\n",sep="")
cmd.out<-cbind(temp1,temp2,temp3,temp4)
temp.cmd<-apply(cmd.out,1,paste,collapse="",sep="")
cat(temp.cmd,file="./cellranger_scripts.txt",append=F,sep="")

There are samples with multiple runs. We run cellranger aggr to merge the multiple runs into one single nUMI vector.

source(file.path(home.dir,"Rprogram/my_functions.R"))
library(xlsx)
home.dir<-"/home/yanxiting/driver_Farnam"
runtable.filepath<-file.path(home.dir,"scratch60/SARC_10X/SraRunTable.txt")
id.filepath<-file.path(home.dir,"scratch60/SARC_10X/scRNA_ids.xlsx")

run.table<-read.table(runtable.filepath,sep="\t",header=T,check.names = F)
id.table<-read.xlsx(id.filepath,sheetIndex = 1,check.names=F)

temp.table<-run.table[run.table$`GEO_Accession (exp)`%in%as.matrix(id.table)[,1],]
temp.list<-split(as.matrix(temp.table)[,1],temp.table$`GEO_Accession (exp)`)

# extract the samples with multiple runs
temp.list<-temp.list[unlist(lapply(temp.list,length))>1]

# find out the fastq part these IDs belong to
temp.filenames<-list.files(file.path(home.dir,"scratch60/SARC_10X"))
temp.filenames<-temp.filenames[grep(".filelist",temp.filenames)]
temp.filenames<-temp.filenames[grep("fastq",temp.filenames)]

file.list<-list()
for(i in 1:length(temp.filenames)){
  temp<-readLines(file.path(home.dir,"scratch60/SARC_10X",temp.filenames[i]))
  temp<-unname(sapply(temp,my.element.extract,splitchar="/",index=-1))
  temp<-unname(sapply(temp,my.element.extract,splitchar="_",index=1))
  temp<-unique(temp)
  file.list[[i]]<-temp
}
names(file.list)<-temp.filenames
temp<-character()
for(i in 1:length(file.list)){
  temp<-c(temp,rep(names(file.list)[i],length(file.list[[i]])))
}

file.list.matrix<-cbind(unlist(file.list),temp)

temp<-list()
for(i in 1:length(temp.list)){
  temp[[i]]<-unique(file.list.matrix[file.list.matrix[,1]%in%temp.list[[i]],2])
}

# transfer the cellranger results back to farnam
# generate the aggregation CSV file
cellranger.dir<-file.path(home.dir,"scratch60/SARC_10X/cellranger_results")
output.dir<-file.path(home.dir,"scratch60/SARC_10X/cellranger_aggr_scripts")

for(i in 1:length(temp.list)){
  output.filepath<-file.path(output.dir,paste(names(temp.list)[i],"_aggr.csv",sep=""))
  # pipeline before cellranger 6.0, use library. Otherwise, use sample_id
  cmd.out<-c("library_id","molecule_h5")
  cmd.out<-rbind(cmd.out,cbind(temp.list[[i]],paste("/home/xy48/scratch60/SARC_10X/cellranger_results/",temp.list[[i]],"/outs/molecule_info.h5",sep="")))
  write.table(cmd.out,file=output.filepath,row.names=F,col.names=F,sep=",",append=F,quote=F)
}

# generate the sh file to run cellranger aggr on the replicated runs of the same sample.
script.dir<-file.path(home.dir,"scratch60/SARC_10X/cellranger_aggr_scripts")
result.dir<-file.path(home.dir,"scratch60/SARC_10X/cellranger_aggr_results")

for(i in 1:length(temp.list)){
  csv.filepath<-file.path(script.dir,paste(names(temp.list)[i],"_aggr.csv",sep=""))
  script.filepath<-file.path(script.dir,paste(names(temp.list)[i],".sh",sep=""))
  cmd.out<-paste("#!/bin/bash\n#SBATCH --time=24:00:00  --ntasks=1 --partition=general --cpus-per-task=1 --mem=80GB --job-name=",names(temp.list)[i]," -o /home/xy48/scratch60/SARC_10X/cellranger_aggr_scripts/",names(temp.list)[i],".sh.o%J -e /home/xy48/scratch60/SARC_10X/cellranger_aggr_scripts/",names(temp.list)[i],".sh.e%J\n",sep="")
  cmd.out<-paste(cmd.out,"module load cellranger/5.0.0\nmodule load bcl2fastq2/2-20-0-foss-2018b\n",sep="")
  cmd.out<-paste(cmd.out,"cd /home/xy48/scratch60/SARC_10X/cellranger_aggr_results\n",sep="")
  cmd.out<-paste(cmd.out,"cellranger aggr --id=",names(temp.list)[i]," --csv=/home/xy48/scratch60/SARC_10X/cellranger_aggr_scripts/",names(temp.list)[i],"_aggr.csv\n",sep="")
  cat(cmd.out,file=script.filepath,append=F)
}

We moved the cellranger aggr results back to cellranger_results together with the samples with unique run.

UMI matrix extraction

We extract the nUMI vector for the samples included in the original paper.

home.dir<-"/home/yanxiting/driver_Farnam"
source(file.path(home.dir,"Rprogram/my_functions.R"))
library(xlsx)
runtable.filepath<-file.path(home.dir,"scratch60/SARC_10X/SraRunTable.txt")
id.filepath<-file.path(home.dir,"scratch60/SARC_10X/scRNA_ids.xlsx")

run.table<-read.table(runtable.filepath,sep="\t",header=T,check.names = F)
id.table<-read.xlsx(id.filepath,sheetIndex = 1,check.names=F)

my.table<-run.table[run.table$`GEO_Accession (exp)`%in%as.matrix(id.table)[,1],]

filenames<-list.files(file.path(home.dir,"scratch60/SARC_10X/cellranger_results"))

merged.data<-numeric()
for(i in 1:length(filenames)){
  if(substr(filenames[i],1,1)=="G"){
    matrix.dir<-file.path(home.dir,"scratch60/SARC_10X/cellranger_results",filenames[i],"outs","count","filtered_feature_bc_matrix")
  }else{
    matrix.dir<-file.path(home.dir,"scratch60/SARC_10X/cellranger_results",filenames[i],"outs","filtered_feature_bc_matrix")
  }
  
  barcode.path <- file.path(matrix.dir, "barcodes.tsv.gz")
  features.path <- file.path(matrix.dir, "features.tsv.gz")
  matrix.path <- file.path(matrix.dir, "matrix.mtx.gz")
  mat <- readMM(file = matrix.path)
  feature.names = read.delim(features.path,
                             header = FALSE,
                             stringsAsFactors = FALSE)
  barcode.names = read.delim(barcode.path,
                             header = FALSE,
                             stringsAsFactors = FALSE)
  if(substr(filenames[i],1,1)=="G"){
    sample.name<-filenames[i]
  }else{
    sample.name<-as.matrix(my.table)[my.table$Run==filenames[i],"Sample Name"]
  }
  
  colnames(mat) = paste0(sample.name,"_",barcode.names$V1)
  rownames(mat) = feature.names$V1
  
  if(i==1){
    merged.data<-mat
    gene.names<-rownames(mat)
  }else{
    merged.data<-cbind(merged.data,mat[gene.names,])
  }
}

output.filepath<-file.path("/home/yanxiting/driver_Farnam/scratch60/SARC_10X/merged_nUMI_45.rds"))

Error: unexpected ')' in "output.filepath<-file.path("/home/yanxiting/driver_Farnam/scratch60/SARC_10X/merged_nUMI_45.rds"))"

We transferred the cellranger results back to farnam to extract the nUMI results. We also add the samples with multiple runs into the merged matrix.

# samples with multiple runs have names of the GEO accession number.
source("/home/yanxiting/driver_Farnam/Rprogram/my_functions.R")
library(Matrix)
cellranger.dir<-"/home/yanxiting/driver_Farnam/scratch60/SARC_10X/cellranger_results"
runtable.filepath<-"/home/yanxiting/driver_Farnam/scratch60/SARC_10X/SraRunTable.txt"

# load in the run table
run.table<-read.table(runtable.filepath,sep="\t",header=T,check.names=F,stringsAsFactors = FALSE)

# load in the run names and the GSM IDs under the cell ranger results folder
temp.filenames<-list.files(cellranger.dir)

sample.rep<-temp.filenames[substr(temp.filenames,1,3)=="GSM"]
sample.uni<-temp.filenames[substr(temp.filenames,1,3)!="GSM"]
my.run.table<-run.table[as.matrix(run.table)[,"Sample Name"]%in%sample.rep | as.matrix(run.table)[,"Run"]%in%sample.uni,]
# extract the nUMI matrix
dirnames<-list.files(cellranger.dir)
dirnames<-dirnames[substr(dirnames,1,4)=="SRR1"]

for(i in 1:length(dirnames)){
  run.name<-dirnames[i]
  sample.name<-run.table[run.table[,"Run"]==run.name,"Sample Name"]
  matrix_dir<-file.path(cellranger.dir,dirnames[i],"outs","filtered_feature_bc_matrix")
  barcode.path <- file.path(matrix_dir, "barcodes.tsv.gz")
  features.path <- file.path(matrix_dir, "features.tsv.gz")
  matrix.path <- file.path(matrix_dir, "matrix.mtx.gz")
  mat <- readMM(file = matrix.path)
  feature.names = read.delim(features.path,
                             header = FALSE,
                             stringsAsFactors = FALSE)
  barcode.names = read.delim(barcode.path,
                             header = FALSE,
                             stringsAsFactors = FALSE)
  colnames(mat) = paste0(sample.name,"_",barcode.names$V1)
  rownames(mat) = feature.names$V1
  
  output.filepath<-file.path("/home/yanxiting/driver_Farnam/scratch60/SARC_10X/cellranger_results",paste0(sample.name,"_",run.name,"_nUMI.rds"))
  saveRDS(mat,file=output.filepath,refhook=NULL)
}

Since we have limited storage on the hpc, we back up the fastq.gz files and the cell ranger output folders. We extracted the nUMI matrix of each sample separately and deleted all other files to make room. All fastq.gz files and their corresponding cell ranger output folders are backed up on google drive under /Volumes/GoogleDrive/My Drive/GraceBackup/GRADS/SARC_PBMC/SARC_10X.

Seurat Analysis

We load the data into Seurat for downstream analysis, especially for the integrated analysis.

First, we create the Seurat object and preliminary filtering on genes and

# due to the unavailability of the cluster, we ran these codes on tac4 server with the same versions of R and packages
# srun --pty --x11 -p pi_kaminski -t 12:00:00 --ntasks=1 --nodes=1 --cpus-per-task=1 --mem=60152 bash
# module restore seurat
#home.dir<-"/home/yanxiting/driver_Farnam"
home.dir<-"/home/yanxiting/Documents/Research/GRADS_SARC_PBMC"
source("/home/yanxiting/Rprogram/my_functions.R")

#output.dir<-file.path(home.dir,"scratch60/SARC_10X/Seurat")
output.dir<-file.path(home.dir,"scRNA-seq/Seurat")

if(file.exists(output.dir)==F){
dir.create(output.dir)
}

# Load the merged data
merged.data<-readRDS(file.path(home.dir,"Data/merged_nUMI_45.rds"),refhook=NULL)

# load in the genes.gtf to convert ensembl ID to gene names
gtf.filepath<-"/home/yanxiting/driver_Farnam/scratch60/SARC_10X/refdata-gex-GRCh38-2020-A/genes/genes.gtf"
gtf.matrix<-read.table(gtf.filepath,sep="\t",header=F,check.names = F, stringsAsFactors = FALSE,skip = 5)
gtf.matrix<-gtf.matrix[gtf.matrix[,3]=="gene",]
temp<-gtf.matrix[,9]

gene.id<-unname(sapply(sapply(temp,my.element.extract,splitchar="; ",index=1),my.element.extract,splitchar=" ",index=-1))
gene.name<-unname(sapply(sapply(temp,my.element.extract,splitchar="; ",index=4),my.element.extract,splitchar=" ",index=-1))
gene.list<-split(gene.id,gene.name)
# get the list of ensembl IDs from mitochondrial genes
mt.genes<-unlist(gene.list[grep("^MT-",names(gene.list))])

#geoff.data <- Read10X(data.dir = file.path(data.dir,"outs","filtered_gene_bc_matrices_mex","GRCh38"))
#load(file.path(home.dir,"scratch_kaminski/public/Backup/Jonas/R_objects/10x_ChuppAsthma.mtx.hybrid_gene_symbols_wo_background_03_1119.Robj"))
gcql1 <- CreateSeuratObject(counts = merged.data, project = "SARC_10X",min.cells = 3,  min.features = 200)
#gcql1[["percent.mt"]] <- PercentageFeatureSet(gcql1, pattern = "^MT-")
gcql1[["percent.mt"]] <- PercentageFeatureSet(gcql1, features = mt.genes)

cat("\n")

cat("Originally, there are ",nrow(merged.data)," genes and ",ncol(merged.data)," cells in the data.\n",sep="")

Originally, there are 36601 genes and 115826 cells in the data.

cat("After the first step filtering, there are ",nrow(gcql1@assays$RNA@data)," genes and ",ncol(gcql1@assays$RNA@data)," cells in the filtered data.\n",sep="")

After the first step filtering, there are 23851 genes and 73145 cells in the filtered data.

rm(merged.data)
#gc(verbose=F)
invisible(gc())

# load in the phenotype data of all the samples
id.filepath<-file.path(home.dir,"Data/scRNA_ids.xlsx")
runtable.filepath<-file.path(home.dir,"Data/SraRunTable.txt")

run.table<-read.table(runtable.filepath,sep="\t",header=T,check.names = F)
id.table<-read.xlsx(id.filepath,sheetIndex = 1,check.names=F)

# get the sample names for all the cells
my.ident<-unname(sapply(colnames(gcql1@assays$RNA@counts),my.element.extract,splitchar="_",index=1))
names(my.ident)<-colnames(gcql1@assays$RNA@counts)

temp<-split(as.matrix(run.table)[,"study_classification"],as.matrix(run.table)[,"Sample Name"])
temp<-lapply(temp,unique)
my.disease<-unlist(temp[my.ident])

my.samplenames<-my.ident
disease.samplenames<-paste(my.disease,"_",my.samplenames,sep="")

# scRNA-seq QC metric.
#mito.genes1 <- grep(pattern = "^MT-", x = rownames(gcql1@assays$RNA@counts), value = TRUE)
#percent.mito1 <- Matrix::colSums(gcql1@assays$RNA@counts[mito.genes1, ])/Matrix::colSums(gcql1@assays$RNA@counts)

# AddMetaData adds columns to object@meta.data, and is a great place to
# stash QC stats
gcql1 <- AddMetaData(object = gcql1, metadata = my.disease, col.name = "disease")
gcql1 <- AddMetaData(object = gcql1, metadata = my.samplenames, col.name = "sample.names")
gcql1 <- AddMetaData(object = gcql1, metadata = my.ident, col.name = "sample")
gcql1 <- AddMetaData(object = gcql1, metadata = disease.samplenames, col.name = "disease.samplenames")
temp<-rep("pbmc",nrow(gcql1@meta.data))
gcql1 <- AddMetaData(object = gcql1, metadata = temp, col.name = "tissue")

# filter the data to only keep genes present (>0 in >=1% of all the cells) in >=3 subjects
gene.pres<-numeric()
sample.names.unique<-unique(gcql1@meta.data$sample.names)

for(i in 1:length(sample.names.unique)){
  temp.matrix<-gcql1@assays$RNA@counts[,gcql1@meta.data$sample.names==sample.names.unique[i]]
  temp.vect<-apply(temp.matrix>0,1,sum)
  gene.pres<-cbind(gene.pres,temp.vect>=(ncol(temp.matrix)*0.01))
}

gene.names<-rownames(gcql1@assays$RNA@counts)[apply(gene.pres,1,sum)>=3]
gcql1.orig<-gcql1
gcql1<-subset(gcql1.orig,features=gene.names)

cat("After the first and second steps of filtering, there are ",nrow(gcql1@assays$RNA@counts)," genes and ",ncol(gcql1@assays$RNA@counts)," cells in the data.\n",sep="")

After the first and second steps of filtering, there are 12834 genes and 73145 cells in the data.

Second, we had a third filtering step on cells based on the nGene and percent.mito. Here are the plots to help us decide the threshold of the filtering.

#g1<-VlnPlot(object = gcql1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,group.by="samplenames",combine=T)

g2<-VlnPlot(object = gcql1, features= c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,group.by = "disease.samplenames",combine=T)
grid.arrange(g2,ncol=1)

# only keep cells with nFeature_RNA>=330 or <=6000
#gcql1<-subset(gcql1,subset= nFeature_RNA >= 300 & nFeature_RNA<= 6000 & percent.mt<= 50)
gcql1<-subset(gcql1,subset= nFeature_RNA >= 300 & percent.mt<= 50)
cat("After the third filtering step, we have ",nrow(gcql1@assays$RNA@counts)," genes and ",ncol(gcql1@assays$RNA@counts)," cells.\n",sep="")

After the third filtering step, we have 12834 genes and 70244 cells.

Based on the quality control plots, in the third filtering step, we removed cells with >50% mitochondrial percentage or <300 expressed genes. After this second filtering step, we have 12834 genes and 70244 cells remained.

Third, we normalize the unimputed and save the SAVER imputed data as the normalized data. These data are saved for downstream further processing.

# normalize the data and save the results
gcql1 <- NormalizeData(gcql1, normalization.method = "LogNormalize",  scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

saveRDS(gcql1, file = file.path(output.dir,"1_normalized_data_unimputed_1.rds"))

#gcql1.saver<-NormalizeData(object = gcql1.saver, normalization.method = "LogNormalize",  scale.factor = 10000)
#saveRDS(gcql1.saver, file = file.path(output.dir,"preprocessed_data_saverimputed_1.rds"))

cat("the preprocessed data for the unimputed data was saved as ",file.path(output.dir,"1_normalized_data_unimputed_1.rds"),"\n",sep="")

the preprocessed data for the unimputed data was saved as /home/yanxiting/Documents/Research/GRADS_SARC_PBMC/scRNA-seq/Seurat/1_normalized_data_unimputed_1.rds

cat("After the normalization, there are ",nrow(gcql1@assays$RNA@data)," genes and ",ncol(gcql1@assays$RNA@data)," cells\n",sep="")

After the normalization, there are 12834 genes and 70244 cells

#cat("the preprocessed data for the SAVER imputed data was saved as ",file.path(output.dir,"preprocessed_data_saverimputed_1.rds"),"\n",sep="")

Original Cell Clustering

First, we find variable genes.

features.n<-2000

We load in the normalized unimputed data and find the top 2000 variable genes for cell clustering. The normalized data was further scaled to remove the effect of nUMI and percent.mito for further cell clustering purpose.

Then we scale the data to remove the effect from nUMI and percent.mito. Here we decide to scale for all genes so that later on, the heatmap will be generated using the scaled data.

#maybe regress cell cycle http://satijalab.org/seurat/cell_cycle_vignette.html
all.genes<-rownames(gcql1)
#gcql1 <- ScaleData(gcql1,features = all.genes, vars.to.regress = c("nCount_RNA", "percent.mt"))
gcql1 <- ScaleData(gcql1)

#output.filepath<-file.path(output.dir,"2_scaled_allgenes_unimputed.rds")
#saveRDS(gcql1, file = output.filepath)

We perform the principal component analysis on the variable genes to decide the number of PCs to use for clustering.

gcql1 <- RunPCA(object = gcql1, features = VariableFeatures(gcql1), npcs=200,verbose=F)

To see if there are PCs that are subject specific, the 2D PCA plots using the top 3 PCs are as follows for the un-imputed data.

#par(mfrow=c(2,2))
fig.1<-DimPlot(object = gcql1, dims = c(1,2), pt.size=0.5)
fig.2<-DimPlot(object = gcql1, dims = c(2,3), pt.size=0.5)
fig.3<-DimPlot(object = gcql1, dims = c(1,3), pt.size=0.5)
grid.arrange(fig.1,fig.2,fig.3,ncol=1)

We draw the heatmap of the top genes for each PC.

DimHeatmap(gcql1, dims = 1:10, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 11:20, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 21:30, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 31:40, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 41:50, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 51:60, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 61:70, cells = 500, balanced = TRUE)

DimHeatmap(gcql1, dims = 71:80, cells = 500, balanced = TRUE)

We use the following Jackstraw plots and the Elbowplot to decide the number of PCs to use for downstream analysis.

gcql1 <- JackStraw(object = gcql1, num.replicate = 100, verbose = TRUE,dims=200)

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~02d 05h 12m 27s
  |+                                                 | 2 % ~02d 05h 54m 46s
  |++                                                | 3 % ~02d 05h 33m 56s
  |++                                                | 4 % ~02d 06h 06m 09s
  |+++                                               | 5 % ~02d 06h 15m 49s
  |+++                                               | 6 % ~02d 05h 50m 36s
  |++++                                              | 7 % ~02d 06h 46m 44s
  |++++                                              | 8 % ~02d 05h 59m 25s
  |+++++                                             | 9 % ~02d 05h 03m 28s
  |+++++                                             | 10% ~02d 04h 46m 39s
  |++++++                                            | 11% ~02d 04h 04m 22s
  |++++++                                            | 12% ~02d 03h 19m 08s
  |+++++++                                           | 13% ~02d 02h 59m 59s
  |+++++++                                           | 14% ~02d 02h 18m 28s
  |++++++++                                          | 15% ~02d 02h 02m 28s
  |++++++++                                          | 16% ~02d 01h 30m 52s
  |+++++++++                                         | 17% ~02d 00h 38m 13s
  |+++++++++                                         | 18% ~02d 00h 03m 54s
  |++++++++++                                        | 19% ~01d 23h 20m 43s
  |++++++++++                                        | 20% ~01d 22h 25m 39s
  |+++++++++++                                       | 21% ~01d 21h 48m 59s
  |+++++++++++                                       | 22% ~01d 21h 11m 54s
  |++++++++++++                                      | 23% ~01d 20h 41m 51s
  |++++++++++++                                      | 24% ~01d 20h 04m 16s
  |+++++++++++++                                     | 25% ~01d 19h 24m 53s
  |+++++++++++++                                     | 26% ~01d 18h 44m 19s
  |++++++++++++++                                    | 27% ~01d 18h 09m 19s
  |++++++++++++++                                    | 28% ~01d 17h 40m 44s
  |+++++++++++++++                                   | 29% ~01d 17h 15m 57s
  |+++++++++++++++                                   | 30% ~01d 16h 45m 35s
  |++++++++++++++++                                  | 31% ~01d 16h 04m 13s
  |++++++++++++++++                                  | 32% ~01d 15h 25m 03s
  |+++++++++++++++++                                 | 33% ~01d 14h 53m 22s
  |+++++++++++++++++                                 | 34% ~01d 14h 18m 09s
  |++++++++++++++++++                                | 35% ~01d 13h 43m 30s
  |++++++++++++++++++                                | 36% ~01d 13h 04m 32s
  |+++++++++++++++++++                               | 37% ~01d 12h 25m 20s
  |+++++++++++++++++++                               | 38% ~01d 11h 59m 17s
  |++++++++++++++++++++                              | 39% ~01d 11h 24m 27s
  |++++++++++++++++++++                              | 40% ~01d 10h 45m 58s
  |+++++++++++++++++++++                             | 41% ~01d 10h 09m 04s
  |+++++++++++++++++++++                             | 42% ~01d 09h 36m 01s
  |++++++++++++++++++++++                            | 43% ~01d 08h 59m 30s
  |++++++++++++++++++++++                            | 44% ~01d 08h 26m 53s
  |+++++++++++++++++++++++                           | 45% ~01d 07h 51m 58s
  |+++++++++++++++++++++++                           | 46% ~01d 07h 13m 28s
  |++++++++++++++++++++++++                          | 47% ~01d 06h 35m 32s
  |++++++++++++++++++++++++                          | 48% ~01d 06h 00m 43s
  |+++++++++++++++++++++++++                         | 49% ~01d 05h 25m 01s
  |+++++++++++++++++++++++++                         | 50% ~01d 04h 54m 35s
  |++++++++++++++++++++++++++                        | 51% ~01d 04h 20m 22s
  |++++++++++++++++++++++++++                        | 52% ~01d 03h 44m 46s
  |+++++++++++++++++++++++++++                       | 53% ~01d 03h 06m 54s
  |+++++++++++++++++++++++++++                       | 54% ~01d 02h 32m 34s
  |++++++++++++++++++++++++++++                      | 55% ~01d 01h 58m 00s
  |++++++++++++++++++++++++++++                      | 56% ~01d 01h 22m 57s
  |+++++++++++++++++++++++++++++                     | 57% ~01d 00h 47m 49s
  |+++++++++++++++++++++++++++++                     | 58% ~01d 00h 11m 25s
  |++++++++++++++++++++++++++++++                    | 59% ~23h 35m 24s  
  |++++++++++++++++++++++++++++++                    | 60% ~23h 01m 55s  
  |+++++++++++++++++++++++++++++++                   | 61% ~22h 28m 22s  
  |+++++++++++++++++++++++++++++++                   | 62% ~21h 54m 21s  
  |++++++++++++++++++++++++++++++++                  | 63% ~21h 20m 02s  
  |++++++++++++++++++++++++++++++++                  | 64% ~20h 45m 15s  
  |+++++++++++++++++++++++++++++++++                 | 65% ~20h 09m 23s  
  |+++++++++++++++++++++++++++++++++                 | 66% ~19h 35m 07s  
  |++++++++++++++++++++++++++++++++++                | 67% ~19h 00m 51s  
  |++++++++++++++++++++++++++++++++++                | 68% ~18h 26m 05s  
  |+++++++++++++++++++++++++++++++++++               | 69% ~17h 50m 29s  
  |+++++++++++++++++++++++++++++++++++               | 70% ~17h 15m 01s  
  |++++++++++++++++++++++++++++++++++++              | 71% ~16h 40m 46s  
  |++++++++++++++++++++++++++++++++++++              | 72% ~16h 05m 28s  
  |+++++++++++++++++++++++++++++++++++++             | 73% ~15h 31m 50s  
  |+++++++++++++++++++++++++++++++++++++             | 74% ~14h 57m 32s  
  |++++++++++++++++++++++++++++++++++++++            | 75% ~14h 23m 28s  
  |++++++++++++++++++++++++++++++++++++++            | 76% ~13h 48m 32s  
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~13h 14m 35s  
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~12h 40m 43s  
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~12h 05m 56s  
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~11h 31m 24s  
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~10h 56m 50s  
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~10h 21m 26s  
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~09h 47m 24s  
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~09h 12m 26s  
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~08h 38m 07s  
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~08h 02m 53s  
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~07h 28m 14s  
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~06h 53m 27s  
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~06h 18m 30s  
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~05h 43m 52s  
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~05h 09m 38s  
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~04h 35m 19s  
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~04h 00m 59s  
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~03h 26m 30s  
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~02h 51m 58s  
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~02h 17m 24s  
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~01h 43m 07s  
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~01h 08m 43s  
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~34m 21s      
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02d 09h 14m 14s

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~06s          
  |+                                                 | 2 % ~05s          
  |++                                                | 3 % ~05s          
  |++                                                | 4 % ~05s          
  |+++                                               | 5 % ~05s          
  |+++                                               | 6 % ~05s          
  |++++                                              | 7 % ~05s          
  |++++                                              | 8 % ~05s          
  |+++++                                             | 9 % ~05s          
  |+++++                                             | 10% ~05s          
  |++++++                                            | 11% ~05s          
  |++++++                                            | 12% ~05s          
  |+++++++                                           | 13% ~05s          
  |+++++++                                           | 14% ~05s          
  |++++++++                                          | 15% ~05s          
  |++++++++                                          | 16% ~04s          
  |+++++++++                                         | 17% ~04s          
  |+++++++++                                         | 18% ~04s          
  |++++++++++                                        | 19% ~04s          
  |++++++++++                                        | 20% ~04s          
  |+++++++++++                                       | 21% ~04s          
  |+++++++++++                                       | 22% ~04s          
  |++++++++++++                                      | 23% ~04s          
  |++++++++++++                                      | 24% ~04s          
  |+++++++++++++                                     | 25% ~04s          
  |+++++++++++++                                     | 26% ~04s          
  |++++++++++++++                                    | 27% ~04s          
  |++++++++++++++                                    | 28% ~04s          
  |+++++++++++++++                                   | 29% ~04s          
  |+++++++++++++++                                   | 30% ~04s          
  |++++++++++++++++                                  | 31% ~04s          
  |++++++++++++++++                                  | 32% ~04s          
  |+++++++++++++++++                                 | 33% ~04s          
  |+++++++++++++++++                                 | 34% ~04s          
  |++++++++++++++++++                                | 35% ~04s          
  |++++++++++++++++++                                | 36% ~04s          
  |+++++++++++++++++++                               | 37% ~03s          
  |+++++++++++++++++++                               | 38% ~03s          
  |++++++++++++++++++++                              | 39% ~03s          
  |++++++++++++++++++++                              | 40% ~03s          
  |+++++++++++++++++++++                             | 41% ~03s          
  |+++++++++++++++++++++                             | 42% ~03s          
  |++++++++++++++++++++++                            | 43% ~03s          
  |++++++++++++++++++++++                            | 44% ~03s          
  |+++++++++++++++++++++++                           | 45% ~03s          
  |+++++++++++++++++++++++                           | 46% ~03s          
  |++++++++++++++++++++++++                          | 47% ~03s          
  |++++++++++++++++++++++++                          | 48% ~03s          
  |+++++++++++++++++++++++++                         | 49% ~03s          
  |+++++++++++++++++++++++++                         | 50% ~03s          
  |++++++++++++++++++++++++++                        | 51% ~03s          
  |++++++++++++++++++++++++++                        | 52% ~03s          
  |+++++++++++++++++++++++++++                       | 53% ~03s          
  |+++++++++++++++++++++++++++                       | 54% ~03s          
  |++++++++++++++++++++++++++++                      | 55% ~02s          
  |++++++++++++++++++++++++++++                      | 56% ~02s          
  |+++++++++++++++++++++++++++++                     | 57% ~02s          
  |+++++++++++++++++++++++++++++                     | 58% ~02s          
  |++++++++++++++++++++++++++++++                    | 59% ~02s          
  |++++++++++++++++++++++++++++++                    | 60% ~02s          
  |+++++++++++++++++++++++++++++++                   | 61% ~02s          
  |+++++++++++++++++++++++++++++++                   | 62% ~02s          
  |++++++++++++++++++++++++++++++++                  | 63% ~02s          
  |++++++++++++++++++++++++++++++++                  | 64% ~02s          
  |+++++++++++++++++++++++++++++++++                 | 65% ~02s          
  |+++++++++++++++++++++++++++++++++                 | 66% ~02s          
  |++++++++++++++++++++++++++++++++++                | 67% ~02s          
  |++++++++++++++++++++++++++++++++++                | 68% ~02s          
  |+++++++++++++++++++++++++++++++++++               | 69% ~02s          
  |+++++++++++++++++++++++++++++++++++               | 70% ~02s          
  |++++++++++++++++++++++++++++++++++++              | 71% ~02s          
  |++++++++++++++++++++++++++++++++++++              | 72% ~02s          
  |+++++++++++++++++++++++++++++++++++++             | 73% ~01s          
  |+++++++++++++++++++++++++++++++++++++             | 74% ~01s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~01s          
  |++++++++++++++++++++++++++++++++++++++            | 76% ~01s          
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~01s          
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~01s          
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~01s          
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05s  

gcql1 <- ScoreJackStraw(gcql1, dims = 1:200)

JackStrawPlot(gcql1, dims = 1:100)

ElbowPlot(object = gcql1,ndims=100)

pc.num<-35

Based on the JackStraw plot and the Elbow plot, we decided to use the top 35 PCs for the unimputed data.

saveRDS(gcql1, file = file.path(output.dir,paste("2_dimreduction_data_",length(gcql1@assays$RNA@var.features),"vargenes_pipeline1.rds",sep="")))
# we save the file as preprocessed_data_saver.rds to prevent any mistakes but we change the file name back to preprocessed_data.rds before performing downstream analysis.

cat("We save the R project as\n")

We save the R project as

cat("unimputed data:\t",file.path(output.dir,paste("2_dimreduction_data_",length(gcql1@assays$RNA@var.features),"vargenes_pipeline1.rds",sep="")),"\n",sep="")

unimputed data: /home/yanxiting/Documents/Research/GRADS_SARC_PBMC/scRNA-seq/Seurat/2_dimreduction_data_2000vargenes_pipeline1.rds

Based on the Jackstraw plots, we decided to use the 35 top PCs for clustering.

Data visualization

UMAP

To identify potential outlying cells and samples, we first label the tSNE plot using disease_samplenames.

sample1<-FindNeighbors(gcql1,dims=1:pc.num)

Computing nearest neighbor graph
Computing SNN

sample1<-FindClusters(sample1,resolution=1.2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 70244
Number of edges: 2352218

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8723
Number of communities: 43
Elapsed time: 24 seconds

1 singletons identified. 42 final clusters.

sample1<-RunUMAP(sample1,dims=1:pc.num)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
10:23:09 UMAP embedding parameters a = 0.9922 b = 1.112
10:23:09 Read 70244 rows and found 35 numeric columns
10:23:09 Using Annoy for neighbor search, n_neighbors = 30
10:23:09 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:23:32 Writing NN index file to temp file /tmp/RtmpqMZeUs/filed37f3287ef
10:23:32 Searching Annoy index using 1 thread, search_k = 3000
10:24:06 Annoy recall = 100%
10:24:08 Commencing smooth kNN distance calibration using 1 thread
10:24:14 Initializing from normalized Laplacian + noise
10:24:21 Commencing optimization for 200 epochs, with 3226940 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:25:55 Optimization finished

fig1<-DimPlot(sample1,reduction="umap")
fig2<-DimPlot(sample1,reduction="umap",group.by="disease.samplenames")
fig3<-DimPlot(sample1,reduction="umap",group.by="disease")
grid.arrange(fig1,fig2,fig3,ncol=1)

tSNE plot

sample1 <- RunTSNE(sample1, dims= 1:pc.num)

Error in Rtsne.default(X = object, dims = dim.embed, ...) : 
  Remove duplicates before running TSNE.

data saving

saveRDS(sample1, file = file.path(output.dir,paste("2_visualization_data_",length(gcql1@assays$RNA@var.features),"vargenes_pipeline1.rds",sep="")))
cat("We saved the seurat object with UMAP and tSNE as ",file.path(output.dir,paste("2_visualization_data_",length(sample1@assays$RNA@var.features),"vargenes_pipeline1.rds",sep="")),"\n",sep="")

We saved the seurat object with UMAP and tSNE as /home/yanxiting/Documents/Research/GRADS_SARC_PBMC/scRNA-seq/Seurat/2_visualization_data_2000vargenes_pipeline1.rds

Integrative Analysis

my.distinct.colors<-c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#ffffff', '#000000')

To describe the cell types captured in the data, the integrated analysis may provide cleaner results. Therefore, we conduct the integrated analysis of all the asthma 10X data using Seurat V3.1 in this note to reduce the subject effect in the data visulization and the trajectory analysis at least.

data loading

########################################################
# 1. load in the singler object
data.dir<-output.dir
sample1<-readRDS(file.path(data.dir,"2_visualization_data_2000vargenes_pipeline1.rds"),refhook = NULL)

sample1.raw<-sample1

integrating data

Split the whole dataset based on subject ID (fresh and DSMO samples from the same subject are considered as one sample). Find anchors to integrate all the datasets together.

sample1.combined<-IntegrateData(anchorset=sample1.anchors,dims=1:40,verbose=FALSE)

Warning: Adding a command log without an assay associated with it

output.filepath<-file.path(output.dir,"integrated_analysis","seurat_object_integrated_allcells.rds")
saveRDS(sample1.combined,file=output.filepath)
cat("We saved the integrated data for all cells as ",output.filepath,"\n",sep="")

We saved the integrated data for all cells as /home/yanxiting/Documents/Research/GRADS_SARC_PBMC/scRNA-seq/Seurat/integrated_analysis/seurat_object_integrated_allcells.rds

output.filepath<-file.path(output.dir,"integrated_analysis","seurat_object_original_allcells.rds")
saveRDS(sample1,file=output.filepath)
cat("We saved the original data for all cells as ",output.filepath,"\n",sep="")

We saved the original data for all cells as /home/yanxiting/Documents/Research/GRADS_SARC_PBMC/scRNA-seq/Seurat/integrated_analysis/seurat_object_original_allcells.rds

Cluster the cells using the integrated data.

DefaultAssay(sample1.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
sample1.combined <- ScaleData(sample1.combined, verbose = FALSE)
sample1.combined <- RunPCA(sample1.combined, npcs = 40, verbose = FALSE)
# t-SNE and Clustering
sample1.combined <- RunUMAP(sample1.combined, reduction = "pca", dims = 1:40)

22:37:42 UMAP embedding parameters a = 0.9922 b = 1.112
22:37:43 Read 69994 rows and found 40 numeric columns
22:37:43 Using Annoy for neighbor search, n_neighbors = 30
22:37:43 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:38:14 Writing NN index file to temp file /tmp/RtmpqMZeUs/filed37f111cb3
22:38:14 Searching Annoy index using 1 thread, search_k = 3000
22:38:53 Annoy recall = 100%
22:38:56 Commencing smooth kNN distance calibration using 1 thread
22:39:02 Initializing from normalized Laplacian + noise
22:39:15 Commencing optimization for 200 epochs, with 3669140 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:41:00 Optimization finished

sample1.combined <- FindNeighbors(sample1.combined, reduction = "pca", dims = 1:40)

Computing nearest neighbor graph
Computing SNN

sample1.combined <- FindClusters(sample1.combined, resolution = 1.2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 69994
Number of edges: 4093949

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8816
Number of communities: 47
Elapsed time: 34 seconds

2 singletons identified. 45 final clusters.

We save the result in a rds file.

output.subdir<-file.path(output.dir,"integrated_analysis")
if(file.exists(output.subdir)==F){
  dir.create(output.subdir)
}

output.filepath<-file.path(output.subdir,"seurat_object_integrated_allcells_cellclustering_noSingleR.rds")
saveRDS(sample1.combined,file=output.filepath)
cat("We saved the seurat object of the integrated data with cell clustering results only as ",output.filepath,"\n",sep="")

We saved the seurat object of the integrated data with cell clustering results only as /home/yanxiting/Documents/Research/GRADS_SARC_PBMC/scRNA-seq/Seurat/integrated_analysis/seurat_object_integrated_allcells_cellclustering_noSingleR.rds

Visualizing integrated data

library(cowplot)
p1 <- DimPlot(sample1.combined, reduction = "umap", group.by = "sample.names")
p2 <- DimPlot(sample1.combined, reduction = "umap", label = TRUE)
#p3<-DimPlot(sample1.combined, reduction = "umap", group.by="singler.hpca.cluster.merged", label = TRUE,cols=my.distinct.colors)
#p4<-DimPlot(sample1.combined, reduction = "umap", group.by="singler.hpca.cluster", label = TRUE,cols=my.distinct.colors)
p5<-DimPlot(sample1.combined, reduction = "umap", group.by="disease", label = TRUE,cols=my.distinct.colors)
#plot_grid(p2, p1, p3,p4,p5,ncol=1)
plot_grid(p2, p1,p5,ncol=1)

Integrate with existing data

We obtained another PBMC 10X data annotated in Kaminski lab so that we can integrate this dataset with the annotated data to annotate the cells in this dataset.

Visualize the integrated data by splitting it to different clusters.

fig.list<-list()
celltype.names<-as.numeric(unique(as.character(Idents(sample1.combined))))
celltype.names<-sort(celltype.names,decreasing=F)
for(i in 1:length(celltype.names)){
  temp.map<-as.character(Idents(sample1.combined))
  temp.map[temp.map!=celltype.names[i]]="other"
  sample1.combined@meta.data$temp.map<-temp.map
  fig.list[[2*(i-1)+1]]<-DimPlot(sample1.combined, reduction = "umap", label = FALSE,group.by="ident")
  #fig.list[[2*i]]<-DimPlot(sample1.combined, reduction = "umap", label = FALSE,group.by="temp.map",cols=c(my.distinct.colors[i],"grey"))
  fig.list[[2*i]]<-DimPlot(sample1.combined, reduction = "umap", label = FALSE,group.by="temp.map",cols=c("red","grey"))
}

do.call(grid.arrange,c(fig.list,ncol=2))

Visualize the integrated data by splitting it to different subjects and different cell type annotation (old).

fig.list<-list()
celltype.names<-unique(as.character(sample1.combined@meta.data$singler.hpca.cluster))
celltype.names<-sort(celltype.names,decreasing=T)
for(i in 1:length(celltype.names)){
  temp.map<-as.character(sample1.combined@meta.data$singler.hpca.cluster)
  temp.map[temp.map!=celltype.names[i]]<-"zother"
  sample1.combined@meta.data$temp.map<-temp.map
  fig.list[[i]]<-DimPlot(sample1.combined, reduction = "umap",group.by="temp.map",label = FALSE,pt.size=0.5,cols=c(my.distinct.colors[i],"grey"))
}
do.call(grid.arrange,c(fig.list,ncol=2))